Hierarchy of N-point functions in the ACDM and ReBEL cosmologies 
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In this work we investigate higher order statistics for the ACDM and ReBEL scalar-interacting 
dark matter models by analyzing 180 /i _1 Mpc dark matter N-body simulation ensembles. The N- 
point correlation functions and the related hierarchical amplitudes, such as skewness and kurtosis, 
are computed using the Count-In-Cells method. Our studies demonstrate that the hierarchical 
amplitudes S n of the scalar-interacting dark matter model significantly deviate from the values 
in the ACDM cosmology on scales comparable and smaller then the screening length r a of a given 
scalar-interacting model. The corresponding additional forces that enhance the total attractive force 
exerted on dark matter particles at galaxy scales lowers the values of the hierarchical amplitudes S„. 
We conclude that hypothetical additional exotic interactions in the dark matter sector should leave 
detectable markers in the higher-order correlation statistics of the density field. We focussed in detail 
on the redshift evolution of the dark matter field's skewness and kurtosis. From this investigation we 
find that the deviations from the canonical ACDM model introduced by the presence of the "fifth" 
force attain a maximum value at redshifts 0.5 < z < 2. We therefore conclude that moderate redshift 
data are better suited for setting observational constraints on the investigated ReBEL models. 

PACS numbers: 98.80.-k, 95.35.+d, 98.65.Dx 



I. INTRODUCTION 

The standard hierarchical structure formation scenario 
assumes that the distribution of mass in the universe has 
grown out of primordial post-inflationary Gaussian den- 
sity and velocity perturbations via gravitational instabil- 
ity. The resulting large-scale structures can be described 
in a statistical way. The two-point and higher order cor- 
relation functions are the most widely studied measures. 
For the standard cold dark matter paradigm - which now 
is a part of the commonly accepted ACDM model - these 
have been studied analytically [e.g. Q]-ii|, as wei l as nu ~ 
merically on the basis of N-body cosmological simulations 

Here we concentrate our study on a modified dark mat- 
ter model that includes long-range scalar interactions be- 
tween DM particles. We focus on the phenomenological 
model of such a long range "fifth" DM force proposed in 
a study by Farrar, Gubser and Peebles |8l-[l3j. We fol- 
low Keselman et at. in dubbing this long-range scalar 
interaction model as ReBEL, daRK Breaking of Equiva- 
lence principLc. This model was proposed as a possible 
remedy for some of the ACDM problems, which relate 
mostly to galaxy scales. For an excellent discussion of 
the motivation behind the long-ran ge s calar-interacting 
model we refer to papers by Peebles |l4 , [ToT | and a recent 
review by Peebles & Nusser [l6||. Over the past few years, 
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the ReBEL model has been extended and explored in a 
range of studies [l7l - [23l |. These studies have revealed its 
potential on the basis of promising results. A variety of 
similar models have also been studied, mostly by means 
of N-body simulations [24l - [3ll ] . There are also additional 
observational arguments in favour of the "fifth force" in 
the dark matter sector, recently forwarded by (32| . 

In this paper we study the hierarchy of N-point cor- 
relation functions of the scalar-interacting DM ReBEL 
model. In principle, these can be used to infer observa- 
tional constraints on the free parameters of the model. 
This is not an entirely trivial affair, since the compar- 
ison of the results with observations is somewhat com- 
plicated by a few factors: (1) galaxies do not necessarily 
trace the mass (biasing) and (2) in the ReBEL model the 
baryonic matter is insensitive to the extra scalar forces. 
Nonetheless, we expect that the information content of 
the higher order correlation functions is sufficient to dis- 
tinguish between the standard DM and scalar-interacting 
DM paradigms. 

To study the high order correlations patterns of the 
DM density field we use cosmological N-body simula- 
tions. The scale and resolution of the simulations are 
designed such that they arc perfectly suited for our pur- 
pose, i.e. they address the highly nonlinear evolution at 
scales smaller than 10 /i _1 Mpc. For the purpose of 
distinguishing between cosmologies these scales are par- 
ticularly useful, since (1) the expected deviation of the 
ReBEL model from the canonical ACDM is maximal at 
these small fully nonlinear scales [Tt], EH, an d (2) nearly 
all detailed observations, except for the largest galaxy 
catalogs, relate to the small or intermediate scales. 

This paper is organized as follows: in section [TT] we de- 
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scribe scalar-interacting DM RcBEL model, followed in 
section Hill by the description of the numerical modelling. 
Section IIVI covers the issues related to the Counts-In- 
Cells method to sample the N-point correlation functions. 
The results of our study are presented in section IVI1 fol- 
lowed by the conclusions in section IVTll 



II. SCALAR-INTERACTING DARK MATTER 
MODEL 

Following our previous work [l8j we study the model 
of the ReBEL long-range scalar interactions in the dark 
matter sector. In this scenario dark matter particles in- 
teract by means of an additional "fifth" force mediated 
by a massless scalar. The extra force term is long-range, 
even though it is dynamically screened by a sea of light 
particles coupled to the scalar field [IS [HI ■ The resulting 
effective gravitational potential between two DM parti- 
cles has the form 113 : 



$ r = 9{x) 

r 

in which G is Newton's constant and 
g{x) = l + pe~ x ^ 



(1) 



(2) 



In this expression r and x are the particle separation in 
real and comoving space. The cosmological scale factor 
a(t) at cosmological time t is normalized to unity at the 
present epoch, a (to) = 1 • The model is specified by 
means of two parameters: 

• /3: strength parameter 

The strength parameter [3 is a dimensionless mea- 
sure of the strength of the scalar interaction with 
respect to a pure Newtonian gravitational gravita- 
tional force: for /3 — 1 the RcBEL forces between 
two dark matter particles are of the same magni- 
tude and strength as the Newtonian gravitational 
force. 

• r s : scale parameter 

the comoving screening length in ft, _1 Mpc, which 
remains constant in the comoving frame. 

The total effective force between two dark matter parti- 
cles of mass mi and ttt-2 is 



F 



DM 



-G 



mi ■ m 2 



1 + 0(1 + —! 



(3) 



From this expression we may immediately infer that the 
regular Newtonian force is recovered at distances r » r s , 
while for separations r < r s the force experienced by the 
dark matter particle will be enhanced or reduced with 
respect to the Newtonian force (depending on the sign of 
the strength parameter 0). 



III. NUMERICAL SIMULATIONS 

A series of N-body numerical experiments is used to 
trace and investigate the growth of the large-scale struc- 
ture in various cosmological scenarios. Part of the simu- 
lations concern the canonical "concordance" ACDM cos- 
mology. Most simulations involve different versions of 
RcBEL cosmologies. In addition, 10 large-scale SCDM 
cosmology simulations are invoked for testing purposes. 

A listing of the parameters and settings of the ensem- 
bles of the simulations is provided by Table HI Simula- 
tions of the concordance ACDM cosmology are labeled 
with LCDM, while the ReBEL ones are labeled with B 
and RS and related parameters indicating the f3 and r s 
parameters of the scalar-interacting dark matter. The 
digits at the beginning of each label relate to the size of 
the simulation box. In addition to the specific scenario 
characteristics - such as fi m , Da, Hubble parameter H, 
us and RcBcl Parameters /3 and r s - the simulations dif- 
fer in terms of the simulation box size Lb ox , number of 
particles N part , force resolution and initial redshift z^. 

The simulations in a 180 ft, _1 Mpc box form the core 
of our study, with the simulations in larger boxes kept 
for additional analysis. With the exception of the 
256LCDMHR and 256B1RS1HR ensembles, all numer- 
ical simulations contain 256 3 dark matter particles to 
sample the theoretical continuum density dark matter 
field. Simulations 256LCDMHR and 256B1RS1HR, con- 
sisting of 512 3 dark matter particles, and simulation en- 
sembles 180LCDMZ80 and 180B1RS1Z80 have a higher 
force resolution, e = 16.8 /i _1 kpc. These simulations are 
used to study the transients and resolution effects. 

For each configuration of simulation parameters we 
generate an ensemble of 5-10 different simulations. This 
enables us to get an estimate of the cosmic variance in- 
troduced by the finite simulation box sizes. Each of the 
ensemble realizations is based on the same amplitude of 
the density field's Fourier components, dictated by the 
power spectrum, while differing in terms of the corre- 
sponding random phases. 

The initial density and velocity fluctuation field in 
all simulations are characterized by a cold dark matter 
spectrum. To generate the initial conditions we use the 
PMcode by Klypin & Holtzman (33[, in conjunction with 
transfer functions computed using the cmbf ast code by 
Seljak & Zaldarriaga [34j| . With the exception of the 
Standard Cold Dark Matter SCDM model, all ACDM 
and ReBEL models start from an initial density field 
with a canonical ACDM power spectrum normalized to 
a linearly extrapolated density variance ag — 0.8 at red- 
shift z = within a sphere of comoving tophat radius 
Rth = 8 h^Mpc. 

The 1024SCDM ensemble traces growth of structure in 
the Standard Cold Dark Matter (SCDM) model. Each of 
the 10 realizations are contained within a 1024 /i _1 Mpc 
cubical box. Even though currently the SCDM model is 
very strongly disfavored by all astronomical data [e.g. [35h- 
l4lj , we use it as reference point and for testing purposes 
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TABLE I: (color on-line) Parameters used in our ensembles of simulations. No. of realizations stands for the number of different 
realizations of the same initial P(k), j3 and r s are scalar-interactions parameters. Lbox denotes the size of the simulation box, 
N P art the number of particles and Zini the initial redshift. The cosmological parameters are: tt m and Qa, denoting the 
dimensionless density parameters of the matter and cosmological constant at redshift z = 0, and as, the amplitude of mass 
fluctuations in a 8 /i _1 Mpc sphere, h is the present dimensionless Hubble parameter, m p is the particle mass, e marks the force 
resolution and I denotes the mean interparticle separation. 
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"These simulations have only 1 realisation, we used 10 bootstrap 
resamplings to obtain the estimates of the mean and standard de- 
viation. 

on the grounds that over thepast decades it has been 
studied in great detail [e.g. EL 0, ■ 

To evolve the particle distribution from the initial scale 
factor to the present time we use the Gadget 2 Tree- 
Particle-Mesh code by Volker Springel (45[, which we 
specifically modified to be able to follow the particle dis- 
tribution in ReBEL force fields. The modifications allow 
the code to handle the long-range scalar-interacting dark 
matter interactions (eqn. [T][3]). The detailed description 
of this modification may be found in our earlier work [l8j . 
Of all simulations, we saved particle positions and veloc- 
ities at redshifts z = 5, 2, 1, 0.5 and 0. The end product 
is a catalog of redshift-dependent snapshots. 

In a simulation with a (comoving) box size of 
180 /i -1 Mpc, a dark matter particle has a mass of 
2.89 x W 10 ^ 1 Mq. In this typical galaxy halo will 

contain roughly a hundred dark matter particles. This 
number is too small to reliably sample any relevant physi- 
cal quantities of a galaxy halo. However, it is sufficient to 
reliably trace the non-linear evolution of the dark matter 
density field down to scales relevant for galaxy formation. 



IV. MOMENTS OF COUNTS-IN-CELLS 

Assuming the applicability of the fair-sample hypothe- 
sis (5fj| , the volume-averaged J-point correlation function 
can be expressed as 

O = Vw [ dx 1 ...dx J iy(x 1 )...iy(x J )0(x 1 ,...,x J ), 

(4) 

where Xi is the comoving separation vector, W(x) is a 
window function with volume 

V w = [ dxH/(x), (5) 
Js 

and the integral covers the entire volume S. Because of 
the fair-sample hypothesis, £j does not depend on the 
location x and is a function of the window volume Vw 
only [H . 



A. Connected Moments 

There is a range of options concerning fast and ac- 
curate methods for measuring the N-point correlation 
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functions of a DM density field sampled by a discrete 
set of N particles. Our analysis is based on the mo- 
ments of the distribution of counts-in-cells (hereafter 

CIC) [HIH El 53- The counts defme a discrete sample 
of the density distribution. Sampling the density field 
by C spherical cells, the J-th central moment of the cell 
counts is defined by 



m 3 (R) 



1 C 



C 



N)- 



(6) 



k 9 = /i 9 - 255fc 2 - 3025k 3 - 7770/c 4 - 6951/€ 5 - 2646/« 6 
- 462/c 7 - 36fc 8 - N. (9) 

Finally the corrected volume-averaged J-th point corre- 
lation functions of DM density field can be written as 



6 = kjN 



-j 



(10) 



; = 1 



We use relations described above to compute £j's up 
to J = 9 from the particle distributions of our N-body 
cosmological simulations. 



where R is the comoving cell radius, iVj the number of 
particles found in a 7-th cell and N the mean number of 
particles in cells of radius R. Following Gaztanagaflrj. 
the connected moments fj,j of the counts may then be 
written as, 



/'2 
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m 2 , 
m 3 , 

o 2 

mi — om 2 , 
m 5 - 10m 3 m 2 , 

TTlQ — 1577147712 



30mS, 



= 7717 — 2177157712 — 3577747773 + 210t773777 2 , 

= 777g — 2877767772 — 567775777,3 — 35777,4 

+ 42O77747772 + 56O77737772 — 630m|, 

= 777g — 3677777772 — 847776?773 — 12677757774 + 75677757772 



252O777477737772 + 560m 



7560777 3 ! 777 2 ; . 



(7) 



The volume-averaged correlation functions £j can be 
computed by dividing the equations for the connected 
moments by N J , 



1; = HjN~ 



(8) 



Shot-Noise effects 



Due to the discrete nature of a finite particle distribu- 
tion, equation [8] is a good estimator of £j only for scales 
where the fluid limit holds. This is satisfied if TV ^ 1. For 
small values of N or, more adequately, for scales compa- 
rable with the mean inter-particle separation, the factor 
/jjN~ J will be dominated by shot noise. 

To correct for the shot noise effects, we use the method 
developed by Gaztanaga [see HB] ■ The method use the 
moment generating function of the Poisson model to cal- 
culate the net contribution by discrete noise. By includ- 
ing this information, one may infer expressions for the 
shot-noise corrected connected moments kj: 



h 

fcg 
hi 

kg 



fMt-N, 

fJ-3 - 3fc 2 - N, 

fJ-4 — 7fc 2 — 6fc 3 - N, 

H 5 - 15fc 2 - 25k 3 - 10fc 4 - N, 

/i 6 - 31fc 2 - 90k 3 - 65fc 4 - 15/l- 5 - N, 

Hr - 63fc 2 - 30U-3 - 350fc 4 - 140fc 5 - 21fc 6 - N, 

/i 8 - 127fc 2 - 966/J3 - 1701fc 4 - 1050/fc 5 - 266k 6 

28fc 7 - N, 



C. Sampling and errors 

Because the computational cost of counting the con- 
tent of cells increases with volume, we adjust the number 
of spherical cells used for the counts-in-cells analysis to 
the comoving cell radius R. We require the total number 
of sampling spheres to be in the range 10 5 < C < 10 6 . 
For the smallest scales we take C = 10 6 , while for the 
largest scales the minimum number of cells is 10 5 . Within 
this range, the number of cells used to sample the mo- 
ments, C(R), scales according to 



C(R) oc 



R 



(11) 



where L is the comoving simulation box width. This 
scaling implies the number of counted points as function 
of scale R to remain comparable. 

Constraining the number of sampling cells is a trade- 
off between the requirement of keeping the sampling er- 
rors as low as possible and limits on the computational 
time. Because the sampling error connected with the fi- 
nite number of cells C scales like C _1 [H]], the decreasing 
number of cells at larger radii R leads to a corresponding 
growth of the intrinsic error. 

In this paper we adopt the standard deviation on the 
mean of the J-point correlation function, determined 
from its estimated values £j in the various realizations 
i (i = 1, . . . , M) within a simulation ensemble (see IJ), 



(£/> 



M 



1 M 



(12) 



as a measure for the variability and error 1 in the esti- 
mate for the correlation function £j, 



Var[6 



\ 



M 



(&-<£/» a , (13) 



The standard deviation of an ensemble obtained by av- 
eraging over its realizations concerns a conservative esti- 
mate of errors. The sampling variance is larger for differ- 
ent realizations within an ensemble than for measurement 
errors associated with the finite number of the sampling 
ceUs El- 
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V. TESTING THE COUNTS-IN-CELLS 
METHOD 

We test our implementation of the CIC method by 
probing its performance with respect to its estimates of 
the two-point correlation function £ 2 and the three-point 
correlation function £3. 



A. Variance and 2nd order moment 

The second order moment is widely used to character- 
ize the rms fluctuation of the matter density field on a 
given scale, 



(14) 



where the scale R is the comoving radius of the applied 
window function W. 

There are two routes towards determining this factor. 
The first estimate of &(R) is yielded by the counts- in- 
cells formalism. Following the Gaztaiiaga formalism, CIC 
leads to the estimate (eqn. QI 



(15) 



£ 2 [CIC](R) ee k 2 N(RY 



where N(R) is the number of particles in spherical cells 
of radius R. 

A second estimate of a(R) is based on the power spec- 
trum P(k) of the dark matter density field in the sim- 
ulations. In theory, the variance follows directly from 
the power spectrum of density fluctuations P(k), via the 
integral over the comoving wave number fc, 



6(i?) = a 2 (R) 



dk 
2^ 



k 2 P(k)W 2 (kR). (16) 



With our analysis being based on counts- in-cells in spher- 
ical volumes of radius R, the natural window function is 
the spherical tophat function. 

In the remainder, the spherical tophat function is used 
as window function. In Fourier space, the top-hat win- 
dow function is specified by 

w run\ o sin ( fcJ ?) - kRcosjkR) 
W TH {kR) = 3 — 3 . (17) 

As a result of the discrete nature of the particles set and 
the finite size of the simulation box, the particle simu- 
lation cannot probe the density perturbations on scales 
larger than the simulation box length L and smaller than 
the mean particle separation, 



I oc N/L 3 



(18) 



(for a simulation of N particles in a box of length L). 
For a proper comparison with the CIC inferred variance, 
the corresponding density field estimate integral in equa- 
tion [TH] is evaluated in between proper integral bound- 
aries. The lower limit is the fundamental mode kr,, while 



the Nyquist frequency k^ yq represents the upper limit. 
For a box of size L, these are 



k L = 



2tt 



ATl/3 



(19) 



where we presume that the number of grid cells on which 
we have sample the initial density field is equal to the 
number of particles N. Hence, the power spectrum vari- 
ance estimate is given by 



i 2 [Pk](R) 



2^ 



2 k 2 P(k)W% H (kR) 



(20) 



For all simulation runs (see table HJ, we have com- 
puted the nonlinear power spectra directly from the 
resulting simulation particle distributions [13]. The 
integral in equation [20] is calculated from the com- 
puted nonlinear power spectra for a limited set of en- 
sembles, those of 1024SCDM, 180LCDM, 180B-05RS1, 
180B02RS1, 180B1RS1 (see tabic ID . 



1. Variance test 

In Fig. [T] we present a comparison between the two 
estimates of the variance a 2 (R), i.e. between the estimate 

£ 2 [CIC] on the basis of the CIC method (eqn. [T5|l and 

the estimate ^[-f^] from the power spectrum integral 
(eqn. [T6|) . For the ensemble of 1024SCDM simulations, 
we determined the variance at three different rcdshifts, 
z = 0,0.5 and z = 1.. The diagram plots the resulting 
variance as a function of the scale R. The symbols (z = 0: 
filled squares, z = 0.5: circles, z = 1: triangles) indicate 
the variance estimates on the basis of the CIC method. 
The continuous lines represent the variance determined 
from the power spectrum integral (z = 0: solid, z = 0.5: 
dashed, z = 1: dotted). 

In the 1024SCDM simulation, the Nyquist frequency 
kNyq ~ 0.785 corresponds to ~ 8 h~ 1 Mpc. This means 
that the diagram in Fig. [1] suffers from a substantial 
level of shotnoise contribution over the range between 
1 /i _1 Mpc < R < 8 /i~ 1 Mpc. Nonetheless, the agreement 
between the two estimators is remarkably good down to 
a scale of w 3 ft, _1 Mpc, comparable to the mean inter 
particle separation in the 1024SCDM ensemble. 



2. Variance Estimate & Model Dependence 

To check whether the modified dynamics of the DM 
fluid in the ReBEL model affects the two variance es- 
timators differently, we compare the resulting estimates 
for a range of different ReBEL models. 

The top panel of fig. [2] compares the two estimates at 
different scales R\y for four different cosmologies: ACDM 
and a ReBEL model with strength parameter (3 = —0.5, 
a ReBEL model with /3 = 0.2 and one with /3 = 1.0. The 
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FIG. 1: (color on-line) Variance Estimators. Comparison of the CIC estimator £ 2 [CIC] and the power spectrum estimator 

£, 2 [Pk] of the variance as a function of scale Rw, based on the ensemble of 1024SCDM simulations. The variance is determined 
for three redshifts: z = 0, z = 0.5 and z — 1.0. The symbols represent the CIC variance estimates. Filled squares: z — 0, 
circles: z — 0.5, triangles: z — 1.. The continous lines indicate the power spectrum integral estimates. Solid line: z = 0, dashed 
line: z — 0.5, dotted line: z = 1. 




FIG. 2: (color on-line) Comparison between the CIC estimator £ 2 [CIC] and the power spectrum integral estimator £ 2 [Pk] of 
the variance of the density field, as a function of scale Rw- The panel shows a comparison between the two estimators for 
simulation ensembles of four different cosmologies. These concern the LCDM cosmology and three different ReBEL cosmologies. 
All simulations have a box size of 180 7i -1 Mpc. The power spectrum estimates are represented by continuous lines, the CIC 
estimates by corresponding symbols. The ensembles are: 1) 180LCDM - LCDM cosmology - solid line - square; 2) 180B-05RS1 
- ReBEL cosmology with P = -0.5 - dotted line - circle; 3) 180B02RS1 - ReBEL cosmology with j3 = 0.2 - dot-dashed line - 
triangle ; 4) 180B1RS1 - ReBEL cosmology - double-dotted line - diamond with /3 = 1.0. For clarity, we only show error bars 
for the 180LCDM ensemble. Top panel: regular plot of variance estimator vs. scale Rw- Bottom panel: Plot of estimator ratio 
i 2 [Pk]/l 2 [CIC]. 
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FIG. 3: (color on-line)The skewness S3 measured for the SCDM simulation ensemble 1024SCDM. The solid line shows the CIC 
estimate of S3 (eqn. [22]) . The dotted line shows the skewness estimated on the basis of perturbation theory (eqn. [23]). The 
error bars (only shown for the CIC) correspond to maximal standard deviation of the ensemble (see main text for the details) . 



lines represent the power spectrum estimate £ 2 Pk ( see 
legend) . The CIC estimates of the variance are indicated 
by symbols of the same colour as the lines, listed in the 
legenda. The difference between the two estimates may 
be best appreciated from the bottom panel, which shows 

the ratio between the two estimators, ^[^^/^[C^C]- 

Both panels clearly shows that for all four different 
cosmologies the two estimators agree very well for scales 
ranging from ps 10 ft, _1 Mpc down to the smallest scales 
that we analyzed, ps 1 /i _1 Mpc. On larger scales, from 
ps 20 ft, _1 Mpc (roughly l/10th of the box width), we see a 
marked disagreement between the two estimators. This 
difference rapidly increases towards larger scales, with 
the CIC estimate systematically increasing as a func- 
tion of scale with respect to the power spectrum value. 
Nonetheless, the fact that the difference between the two 
estimates is identical for the different model ensembles, 
in terms of character and scale at which they start to 
diverge, indicates that the accuracy achieved by the CIC 
method is the same for each of the cosmologies. 



B. Third order moment: the S3 test 

The second order density field statistic, represented 
by the two-point correlation function, is not sufficient for 
characterizing the density field beyond the linear phase of 
structure evolution. Moving into the quasi-linear phase, 
we start to discern the gravitational contraction of over- 



dense regions into sheetlike and filamentary patterns and 
compact dense haloes and the volume expansion of low 
density void regions. To be able to follow and charac- 
terize this process, we need to turn to the higher order 
moments of the density field. 

To test the performance of the CIC estimator, we turn 
to the reduced third moment of the density field. The 
skewness S3 is defined as 



4 



(21) 



An estimate of S3 can therefore be readily obtained 
on the basis of the corrected volume-averaged 3-point 
correlation function of the dark matter density field (see 
eqn. [TU] and eqn. [5]) , 



S 3 [CIC\(R) 
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a 



k 3 N(R)- 
k 2 2 N(RY 



= ^N(R) 7 (22) 



where N(R) is the number of particles in spherical cells 
of radius R. 

An alternative estimate of the skewness finds its origin 
in weakly nonlinear perturbation theory (PT, [l|, 0, S3, 
Hil). Juszkiewicz et a/.Q showed that a good approxi- 
mation for the skewness S3 of the field, smoothed with 
the spherical top-hat window, is given by 
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FIG. 4: (color on-line) Comparison of the measured skewness Sa(R) in ensembles with different box sizes. Top panel: results for 
pure ACDM model simulations. Bottom panel: results for the ensembles for a ReBEL model with scalar interaction parameters 
j3 = 1 and r s = 1 /i _1 Mpc. In each panel we plot the S3 (71) relation for three different ensembles of the same model, one in 
a box with a width of 180 /i _1 Mpc (blue dotted line), of 360 h~ 1 Mpc (green dashed line) and one of 512 ft Mpc (red solid 
line). The vertical lines show the corresponding values of the Nyquist scale Rxyq = 2n/kN yq ~ 21 for each of these simulation 
boxes. The error bars represent the la errors in the 180LCDM (top panel) and 180B1RS1 (bottom panel) ensembles. 



where 71 is the logarithmic slope of the variance, defined 
as 

d log o 2 (R) 

7i = ~(" + 3)= dlQgR , (24) 

where n is the slope of the power spectrum at scale R. 
The term 34/7 is a well-known result pertaining to the 
unsmoothed field (see [l|). For the estimate of the skew- 
ness Ss[PT] based on this result, we use the estimate of 
the variance a 2 (R) obtained via the integral over the non- 
linear power spectrum P(k) for a tophat filter W(kR), 

ie. from %[Pk] (R) (eqn.HJJ), 

§ s[ PT] W - f + " l0 jfJ iR) . (25) 

In fig. [3] we have compared the two estimates of the 
skewness S 3 for the SCDM simulations in the 1024SCDM 
ensemble, over a range of 1 ft _1 Mpc < R < 80 h~ 1 Mpc. 
The solid line represents the skewness measured directly 
on the basis of the counts in cells method (eqn. |2"2"|) . 
while the dotted line is the perturbation theory predic- 
tion (eqn. [33]). The error bars, here shown for the CIC es- 
timates, are the maximal standard deviation of the mea- 
surements for the simulations in the 1024SCDM ensemble 



(taking into account that at each different scale R we use 
a different number of sampling cells). 

Overall, we find that the two skewness estimators are 
in reasonable agreement with each other, in particular 
on linear and mildly nonlinear scales, 8 ft~ x Mpc < R < 
80 ft _1 Mpc, exactly as expected and reported by many 
other authors [I HSU. 



1. The box size test 

The effects of finite volume on the statistics of large 
scale structure have been extensively studied in several 
studies [5fj| . For most of the results presented in this 
paper, finite volume effects arc rather unimportant. We 
focus mainly on a direct comparison between observables 
of the canonical ACDM model and those of the scalar- 
interacting dark matter ReBEL models. As long as any 
of the finite volume induced artefacts affects each of the 
cosmological models to a comparable extent, we need not 
worry about their influence on the results of our study. 

Nonetheless, there is one factor which needs to be in- 
vestigated in some detail. The new physics of the dark 
sector scalar-interacting ReBEL models involves a new 
fundamental and intrinsic scale, the screening length r s . 
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It is a priori unclear in how far the relation between the 
length L of the simulation box and the screening length 
r s of the RcBEL model will be of influence on the counts- 
in-cell measurement of various moments. 

To evaluate whether the finite box size has any im- 
pact on the measured values of 53, we have run a set of 
simulations for two different models. One model is the 
ACDM model, whose gravitational force law is entirely 
scale-free, while the other model is a RcBEL model char- 
acterized by an intrinsic force scale. We chose a RcBEL 
model with strength factor (3 = 1 and scale parame- 
ter r s = 1 /i _1 Mpc. Each of the two sets of simula- 
tions contain three ensembles of the same cosmological 
model. The first ensemble has a 180 /i _1 Mpc simulation 
box, the second a 360 ft. _1 Mpc box and the third one a 
512 h^Mpc box. 

In figure [4] we follow the trend of the skewness Ss(R) 
as a function of scale R, for each of the simulation ensem- 
bles. The top panel shows the results for the three sets 
of ACDM simulations, the 180LCDM simulations in a 
180 fr, _1 Mpc box (dotted line), the 360LCDM simulations 
in a 360 /i _1 Mpc box (dashed line) and the 512LCDM 
simulations in a 512 /i -1 Mpc box (solid line). The same is 
repeated for the ReBEL model in the bottom panel, with 
the 180B1RS1 simulations in a 180 /i _1 Mpc box (dot- 
ted line), the 360B1RS1 simulations in a 360 /i~ 1 Mpc 
box (dashed line) and the 512LCDM simulations in a 
512 /i _1 Mpc box (solid line). In the figure we have 
also indicated the location of the Nyquist scale, RNyq = 
2n /kptyq, of each of the three simulation boxes. The three 
vertical lines mark their position. 

We find that in the ACDM case, the measured skew- 
ness in the simulation ensembles with different box size 
agree very well over the entire ranged we probed, from 
the largest measured scales ~ 30 /i _1 Mpc < R < 
80 /i _1 Mpc down to the smallest scales of 1 /i _1 Mpc < 
R < 8 ft,~ 1 Mpc. Interestingly, we also find a similar 
good agreement between the simulation ensembles of the 
RcBEL model. Moreover, we also find a surprisingly 
good agreement at scales where we expect two-body ef- 
fects to start to dominate, below the Nyquist scale of the 
simulation. 

Given the fact that the measured S3 values remain 
consistent over such a wide range of scales and seems 
independent of the size of the simulation box size, we 
conclude that the effect of a different ratio r s /L of in- 
trinsic force scale to box size has negligible, if any, effect 
on the measurement of statistical moments. 



2. Transients 

The PMcode that we use to generate the initial con- 
ditions is based on the Zeldovich Approximation (ZA) 
method [5l|. It is well known that the Zeldovich approx- 
imation introduces an artificial level of skewness and ad- 
ditional higher order hierarchy moments into the density 
field (52l . [53}. A sufficient number of simulation time- 



steps is required for the true particle dynamics to take 
over and to relax these transient artifacts. An alternative 
approach is to resort to second order Lagrangian pertur- 
bation theory schemes for setting up the initial conditions 
of simulations [52T - [55| . 

Because of the above, the initial redshift of a cosmolog- 
ical simulation is an important factor in determining the 
statistical reliability of the cosmological numerical exper- 
iment. In general, for the purpose of comparing density 
fields and cumulants in different models we need to be 
less concerned about the net amplitude of the transients 
as they will have the same magnitude in all models. 

Nonetheless, there is an additional factor that depends 
on the initial redshift and which only affects the ReBEL 
models. The intrinsic scalar force of these models should 
be able to act as long as possible, in order to account 
for an optimal representation of their impact on the 
dark matter density field. If the ReBEL simulations are 
evolved too far by means of the Zeldovich approximation 
and their dynamical evolution started too late, the devi- 
ation of the ReBEL dark matter density field from the 
one in the ACDM simulations will diminish. 
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FIG. 5: (color on-line) Test of transients effects for 
180 7i _1 Mpc ensembles. S 3 (zi = 40)/S 3 (zj = 80): ratio 
of skewness S3(zi = 40) measured for the simulation en- 
semble started at redshift Zi = 40 to the skewness of the 
simulations started at redshift Zi = 80. The ratio S3 (Zi = 
AQ)/Si{zi — 80) is plotted as a function of scale R. The 
solid black horizontal line represents the unity line for which 
Sz(zi = 40) = Ss(zi = 80). The dashed vertical line marks 
the Nyquist scale for the simulations in a 180 /i _1 Mpc box, at 
27v/k Nyq = 1.38 /i _1 Mpc. The ratio S 3 (zi = 40)/S , 3 (z 4 = 80) 
is plotted for two different situations. Red dashed line: 
skewness ratio for the two ensembles of ACDM simulations, 
180LCDM and 180LCDMZ80. Blue dotted line: skewness ra- 
tio for the two ReBEL simulation ensembles with (3 = 1 and 
r s = 1 ft^Mpc, 180B1RS1 and 180B1RS1Z80. The error 
bars represent the la values determined for the 180LCDM 
ensemble. The errors in the other ensembles have a similar 
magnitude. 
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In order to quantify the possible effects of the tran- 
sients, we have performed a series of auxiliary simulation 
ensembles. These contain 256 3 DM particles placed in 
boxes of the box width 180 /i _1 Mpc and have 10 times 
better force resolution. There are two ensembles, one for 
the ACDM model, 180LCDMZ80, and one for the ReBEL 
model with /3 = 1 and r s = 1 h^Mpc, 180B1RS1Z80. 
We will compare them with our main ensembles for the 
same models, 180LCDM and 180B1RS1. Therefore the 
ensembles of each model will differ only in the force res- 
olution and the redshift at which the N-body calculation 
is started, one at z, ; = 80 and the other at Zi = 40. 

The results for the direct comparison of the skewness 
S3 in the Zj = 80 models and the z,; = 40 models, in 



terms of their ratio 5| i= °/5| i= , are plotted in hgureJS] 
The blue dotted line represents the ratio for the ReBEL 
ensemble, the red dashed line for the ACDM ensemble 
(for which (3 = 0). For reference, the black horizontal 
solid line indicates the unity ratio S^ i=4 ° / S^ i=80 = 1, 
while the vertical line marks the Nyquist scale 2n/kff yq = 
1.4 ft, _1 Mpc for these simulations. The error bars mark 
la errors in the 180LCDM ensemble, with errors in the 
other three ensembles being of the same order. 

We note that the visible transients effects arc, if real, 
very small. The skewness ratio curves lie very close to 



Zi=80 



1. Their deviations from 



the unity line 5| i=40 /5; 
unity are smaller than the la errors, with discrepancies 
not exceeding the 10% level. On the basis of this we may 
conclude that the redshift of the initial conditions of our 
main ensembles, at Zi =40, is sufficiently high to assure 
that any effects of possible transient are negligible for our 
analysis. 



3. Resolution 

The last important effect we must check is the im- 
pact of the mass and force resolution used in our simula- 
tions on the measured quantities. The mass resolution is 
related to the mean inter-particle separation, while the 
force resolution corresponds to the scale at which the 
force prescription of the simulation code exactly recovers 
the intended Newtonian - or ReBEL - force. 

To investigate the impact of these resolution factors 
on the measurement skewness we use the high force res- 
olution ensembles 180LCDMZ80 and 180B1RS1Z80, as 
well as two single high mass and force resolution runs, 
256LCDHR and 256B1RS1HR (see table 0}. For these 
two simulations we use bootstrap resampling to obtain 
averages of mean and variance of the measured moments. 
This is accomplished as follows. We randomly cast a large 
number of spherical cells over the entire simulation vol- 
ume. This ranges from 2 x 10 s cells with R = 1 h Mpc 
to 2 x 10 6 cells for R = 30 h -1 Mpc. Ten sets of mea- 
surements were constructed, each consisting of a random 
subset of 10% of the casted spheres. 

We may assess the resolution effects on the basis of the 
plot in figure [5] It depicts the ratio of the skewness in 
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FIG. 6: (color-on-line) Test of resolution effects on computed 
skewness S3. Plotted are the ratios Si 1RS1 /S£ CDM of the 
skewness obtained in the ReBEL model with )3 — 1 and 
r s = 1 h~ Mpc to the skewness measured in the LCDM 
model, for three simulations with different force and mass 
resolutions. The lines depict relevant ratios for lower reso- 
lution run 180B1RS1 (solid line), high force resolution run 
180B1RS1Z80 (dashed line) and high mass and force resolu- 
tion run 256B1RS1HR (dotted line). The vertical dashed line 
marks Nyquist scale for runs with 256 particles, while the 
vertical dotted line depicts the same scale for 512 3 particles 
simulations. The error bars represent the la values deter- 
mined for the 180B1RS1 ensemble. 



three different ReBEL ensembles to that of the skewness 
in the LCDM model, Si> 1RS1 /S£ CDM . Each of the three 
ReBEL models have the same ReBEL parameters, (3 = 
1 and r s = 1 ft, _1 Mpc, but differ in resolution. The 
lower resolution run is 180B1RS1 (solid line), the high 
force resolution run is 180B1RS1Z80 (dashed line), while 
the high mass plus high force resolution run is that of 
256B1RS1HR (dotted line). The error-bars marking the 
skewness ratio of the 180B1RS1 run are the la errors for 
180B1RS1 ensemble. 

Even though we find that the simulations with a higher 
resolution show a systematically higher signal level at 
scales R < 7 /i _1 Mpc, this effect is entirely contained 
within - or at best marginally above - the la errors of the 
180LCDM ensemble. We may therefore conclude that 
an increase in the force and/or mass resolution of the 
simulation does not yield a significant improvement of 
the signal level. This reassures us that the simulation 
ensembles used in our main study yield good and reliable 
estimates of the quantities which we study. 



C. CIC test summary 

In all, we may conclude from the various tests of the 
Count-in-Ccll method that it is perfectly suited for study- 
ing the impact impact of long-range scalar interactions on 
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the higher-order correlation statistics of the dark matter 
density field. 



VI. MOMENT ANALYSIS OF N-BODY 
ENSEMBLES 

Having ascertained ourselves of the reliability of the 
CIC machinery, we will present and discuss the results of 
the correlation function analysis of our N-body experi- 
ments. The intention of this study is the identification of 
discriminative differences between the canonical ACDM 
cosmology and a range of scalar interaction ReBEL mod- 
els. 

We address two aspects of the resulting dark matter 
distributions. The first concerns a complete census of the 
hierarchy amplitudes S n , from n = 3 to n — 8 for a set of 
three different ReBEL model simulations and a similar 
ensemble of ACDM simulations. In addition, in order to 
assess the redshift evolution of these statistical measures, 
we focus on the redshift dependence of the skewness and 
kurtosis. 



A. Hierarchy amplitudes 

The hierarchy amplitudes S n (R) of order n are con- 
ventionally defined as, 

S n (R) = -^T = £ n *- 2(n - 1) , (26) 
£2 

with the volume-averaged correlation functions £ n (R) 
and variance a 2 (R) implicitly depending on the scale R. 



1. General Trends 

In figure [7] and [5] we plot the measured S„s, from n = 3 
up to n = 8, for all simulation ensembles with boxwidth 
180 /i _1 Mpc (see table HJ. These two figures represent 
the key result of this study. 

The volume-averaged N-point correlation functions £ n 
have been computed by means of the CIC method, fol- 
lowing the description in section IIVI The simulations 
for which the hierarchy amplitudes have been computed 
are the 180LCDM set of ACDM simulations and the 
180B-05RS1, 180B02RS1 and 180B1RS1 simulations of 
the ReBEL models with scalar interaction scale param- 
eter r s = 1 /i -1 Mpc and strength parameter j3 = —0.5, 
f3 = 0.2 and /? = 1.0. The (3 = —0.5 case, whose physical 
effect is that of a repulsive scalar ReBEL force, does not 
have a real physical motivation. It is mainly included for 
reference, in order to outline the impact of the /3 strength 
parameter on the final nonlinear density field. For all 
model simulations we have calculated the hierarchy am- 
plitudes S n at 12 logarithmically spaced scales within 



the range of 1 /i -1 Mpc < R < 36.09 /i _1 Mpc. The exact 
values of these 12 scale values R are listed in table UT1 

Figures [7J and [8] plot the hierarchy amplitudes S n as a 
function of scale R. The two figures are complementary: 
in figure[H]the S„ are shown separately for each of the cos- 
mological models, while figure [7J superimposes the curves 
for each of the models in order to highlight their differ- 
ences. In addition, to provide an impression of the rela- 
tive differences between hierarchy amplitudes in each of 
the cosmological models, figure |H] plots the ratio between 
the S n (R) between each ReBEL model and the concor- 
dance ACDM models. In figure [7J each cosmological 
model is indicated by a different line types. The canon- 
ical ACDM model is indicated by the black solid line, 
the ReBEL model with /3 = -0.5 and r s = 1 /i _1 Mpc 
by the blue dotted line, the ReBEL model with j3 = 0.2 
and r s = 1 /i _1 Mpc by the green dot-dashed line and 
the ReBEL model with f3 = 1.0, r s = 1 fc. _1 Mpc by 
the red dashed line. Figure [8] also includes the error 
bars of the measured S n values, restricted to their up- 
per half for purposes of clarity. For reference, we have 
listed the values of the standard deviation for the skew- 
ness S3 and kurtosis 54 in tables [TT] and IIIII The thin 
dashed vertical lines in figures [7J and [5] mark the Nyquist 
scale RNyq = 2n/k Nyq ps 1.4 /i _1 Mpc for the ACDM and 
ReBEL simulations, which for these realizations is dou- 
ble the mean inter-particle separation 21. We consider the 
computed quantities on scales below the Nyquist scale as 
unreliable, and exclude them from further analysis in this 
study. 

There are some clear trends in the behaviour of the S n 
hierarchy. At large scales, R > 10 ft. Mpc, all cosmolo- 
gies agree on the S n . This is straightforward to under- 
stand because at these large scales the ReBEL models 
are practically equivalent to the ACDM cosmology. The 
differences between the models become distinct at scales 
R < 10 /i _1 Mpc, where the effect of the scalar ReBEL 
force kicks in. We discern a systematic trend, with all S n 
consistently higher than the ACDM values for the ReBEL 
model with j3 = —0.5, consistently lower than the ACDM 
values for the ReBEL model with (3 = 1.0 and the val- 
ues for the ReBEL model with /3 = 0.2 straddling tightly 
around the ACDM values. We also notice that the dif- 
ferences between the models increase systematically as 
a function of order n (see fig. This may be easily 
understood from the higher sensitivity of the higher mo- 
ments to the changing shape of the density probability 
function, and hence to the changes in the dark matter 
density distribution. 



2. Skewness and Kurtosis 

In the observational reality, beset by various sources 
of noise, it may be cumbersome to get reliable estimates 
of higher order moments. On the other hand, we may 
expect reasonably accurate estimates of the third and 
fourth order moments, the skewness and kurtosis. The 
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FIG. 7: (color on-line) Hierarchical amplitude S n (R) as a function of scale R. Plotted are S n , for n = 3 to n = 8, for 
four different simulation ensembles: the canonical ACDM model 180LCDM simulations (solid black line), the /3 = —0.5 and 
r s = 1 /i _1 Mpc ReBEL model simulations 180B-05RS1 (blue dotted line), the P = 0.2 and r s = 1 /i -1 Mpc ReBEL model 
simulations 180B02RS1 (green dot-dashed line) and the fj = 1.0 and r s — 1 ft _1 Mpc ReBEL model simulations 180B1RS1 (red 
dashed line). The 5*3 curves have the lowest amplitude, with the amplitude of the S n curves systematically increasing as a 
function of n. The thin vertical dashed line marks the Nyquist scale RN yq ~ 1.4 ft _1 Mpc. 




FIG. 8: (color on-line)Hierarchical amplitudes S n (R) as a function of scale R. Equivalent to fig. [7] the S n are plotted in a 
separate panel for each cosmological model. Top left: f3 — —0.5 and r s — 1 ft~ 1 Mpc ReBEL model simulations 180B-05RS1; top 
right: (3 = 1.0 and r s = 1 h -1 Mpc ReBEL model simulations 180B1RS1; bottom left: the canonical ACDM model 180LCDM 
simulations; bottom right: /3 = 0.2 and r s = 1 h, _1 Mpc ReBEL model simulations 180B02RS1. The S3 curves have the lowest 
amplitude, with the amplitude of the S n curves systematically increasing as a function of n. The thin vertical dashed line 
marks the Nyquist scale 7?jv a9 ~ 1.4 ft _1 Mpc. 
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FIG. 9: (color on-line) Hierarchical amplitude S n {R) as a function of scale 7?: ratio S n (R) /S n (7?) of the hierarchy amplitudes 
in any of the simulation ensembles to the hierarchy amplitudes in the concordance ACDM cosmology simulations. From top to 
bottom panel: S 3 /5 , 3 LCDM to S , 8 /5'g LCDM . In each panel we plot the curves for four different simulations ensembles: the canonical 
ACDM model 180LCDM simulations, by definition equal to unity (solid black line), the /3 = —0.5 and r s = 1 h~ J Mpc ReBEL 
model simulations 180B-05RS1 (blue dotted line), the /3 = 0.2 and r s = 1 /i _1 Mpc ReBEL model simulations 180B02RS1 
(green dot-dashed line) and the /3 = 1.0 and r s = 1 ft _1 Mpc ReBEL model simulations 180B1RS1 (red dashed line). The S3 
curves have the lowest amplitude, with the amplitude of the S n curves systematically increasing as a function of n. The thin 
vertical dashed line marks the Nyquist scale RN yq ~ 1-4 /i _1 Mpc. The error-bars correspond to la scatter of the 180LCDM 
ensemble. 



14 



TABLE II: Measured values of the S3 hierarchy amplitude 
at redshift z = 0, for four simulation ensembles in differ- 
ent cosmologies. The models are the ACDMmodei 180LCDM 
simulation, the /3 = —0.5 and r s = 1 ft~ 1 Mpc ReBEL model 
simulations 180B-05RS1, the /3 = 0.2 and r s = 1 ft^Mpc 
ReBEL model simulations 180B02RS1 and the fi = 1.0 and 
r s = 1 /i^Mpc ReBEL model simulations 180B1RS1. 



R 180LCDM 
( h^Mpc) 



180B-05RS1 180B02RS1 180B1RS1 



01.00 
01.38 
01.92 
02.66 
03.68 
05.10 
07.07 
09.80 
13.57 
18.80 
26.05 
36.09 



6.53 ±0.35 
6.82 ± 1.23 
6.10 ±0.69 
5.70 ±0.66 
5.05 ±0.85 
4.45 ± 0.46 
3.97 ±0.50 
3.52 ±0.46 
3.15 ±0.48 
2.89 ±0.65 
2.66 ±0.86 
2.60 ± 1.10 



7.65 ± 1.08 
7.69 ± 1.23 
7.29 ± 1.45 
6.35 ±1.37 
5.48 ± 0.87 
4.72 ±0.81 
4.07 ±0.61 
3.56 ±0.50 
3.21 ±0.54 
2.90 ±0.66 
2.68 ±0.87 
2.60 ± 1.10 



6.40 ± 1.00 
6.35 ± 0.74 
6.09 ± 0.70 
5.66 ±1.00 
4.80 ± 0.39 

4.49 ± 0.73 
3.90 ±0.51 

3.50 ± 0.47 
3.12 ± 0.47 
2.87 ±0.65 
2.66 ±0.87 
2.57 ± 1.08 



5.30 ± 0.59 
5.13 ±0.47 
5.35 ±0.84 
4.90 ±0.52 
4.51 ±0.55 
4.26 ±0.51 
3.74 ±0.40 
3.41 ±0.43 
3.10 ±0.47 
2.85 ±0.63 
2.63 ±0.84 
2.60 ±1.09 



TABLE III: (color on-line) Measured values of the £4 hierar- 
chy amplitude at redshift z = 0, for four simulation ensembles 
in different cosmologies. The models are the ACDMmodel 
180LCDM simulation, the /3 = -0.5 and r s = 1 h^Mpc 
ReBEL model simulations 180B-05RS1, the /3 = 0.2 and 
r s = 1 h^Mpc ReBEL model simulations 180B02RS1 and 
the /3 — 1.0 and r s — 1 /i~ 1 Mpc ReBEL model simulations 
180B1RS1. 

R 180LCDM 180B-05RS1 180B02RS1 180B1RS1 
( h^Mpc) 



Assessing the data presented in these tables reveals 
that S3 and S4 values converge to within ler around the 
ACDM values for scales larger than R = 9.8 h^ 1 Mpc. 
As we turn towards smaller scales R, the ReBEL model 
values for the skewness and kurtosis display an increas- 
ingly large difference with respect to the ACDM value. 
In other words, at these small (mildly) nonlinear scales 
we observe a direct imprint of the scalar forces on the 
density field moments. 

At scales comparable to the screening length, ReBEL 
models with a positive strength parameter (3 have a lower 
skewness and kurtosis value than those for the canoni- 
cal ACDM model. The difference is smaller for ReBEL 
models with a lower f3, and turns into a higher value as f3 
turns negative. Seen as a function of scale, the difference 
decreases towards larger scales R. 

For the f3 = 1.0 180B1RS1 simulations the value of S 3 
at R ~ 1.4 h^ 1 Mpc is f=a 25% lower than the value for 
the ACDM model, while the discrepancy is in only the 
order of « 10% at R = 3.68 h" 1 Mpc and has dropped 
towards < 5% for R > 7 h~ 1 Mpc. The differences are 
more prominent in the case of the kurtosis S4. For R ~ 
1.4 ft^Mpc the value of S 4 is smaller than the ACDM 
value by no less than « 48%, decreasing towards « 27% 
at R = 3.68 and to less than 10% at R > 7 /i^Mpc. 

The differences between the cosmological models are 
therefore less substantial for the skewness S3 than for the 
kurtosis S4. On the condition that it is possible to obtain 
reliable estimates for S4 in the observational reality, this 
leads us to the conclusion that the kurtosis may be better 
suited as tracer of ReBEL signatures in the density field. 
More detailed studies and simulations, including baryons 
and mock galaxy samples, will be necessary to make a 
final choice for the optimal marker of ReBEL cosmology 
in observational catalogues. 



Redshift evolution 
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8±5 



question is whether the presence or absence of ReBEL 
scalar forces may be deduced from the behaviour of these 
moments. To evaluate the discriminatory powers of S3 
and S4 we list the measured values of these hierarchy 
amplitudes in tables HTl and IIIII 



In our previous study |18[ we found that the amplitude 
of the deviation of the two-point correlation function £ of 
the ReBEL model to that of the ACDM model changes 
with redshift. This suggests a similar evolution of higher 
order moments like S3 and S4, prodding us to assess the 
redshift evolution of skewness and kurtosis. 

To this end, we study the archive of five snapshots - 
at redshifts z = 5., 2., 1., 0.5 and z = - which for each 
simulation in the four 180 h~ 1 Mpc ensembles were saved: 
180LCDM, 180B-05RS1, 180B02RS1, and 180B1RS1. 

The redshift evolution of the skewness and kurto- 
sis in the four simulation ensembles can be followed 
in figure QT] In the lefthand column we plot the ra- 
tio S 3 HeBEL /S 3 ^ CDM of the skewness in the three differ- 
ent ReBEL models to the one for the canonical ACDM 
model in a sequence of five panels, for the five subse- 
quent redshift snapshots which we analyzed, from z = 5. 
(bottom) to z = (top). The righthand column is 
organized in an equivalent manner for the kurtosis ra- 
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FIG. 10: (color on-line) Skewness Ss(R) and kurtosis S<i(R) 
as a function of scale R. Plotted are Ss(R) and S4,(R) in four 
180 /i _1 Mpc box simulation ensembles: the canonical ACDM 
model 180LCDM simulations (solid black line), the f3 — —0.5 
and r B = 1 /i -1 Mpc ReBEL model simulations 180B-05RS1 
(blue dotted line), the /3 = 0.2 and r s = 1 h _1 Mpc ReBEL 
model simulations 180B02RS1 (green dot-dashed line) and 
the ft — 1.0 and r B = 1 ft - Mpc ReBEL model simulations 
180B1RS1 (red dashed line). The thin vertical dashed line 
marks the Nyquist scale RN yq ~ 1.4 7i _1 Mpc. The error-bars 
correspond to the la scatter of the 180LCDM ensemble. 



tio S? eBEL /S£ CUM . In the panels we follow the same 
nomenclature and line scheme as in the previous sec- 
tion^): the canonical 180LCDM ACDM model is indi- 
cated by the black solid line, the 180B-05RS1 ReBEL 
model with /3 = —0.5 and r s = 1 /i -1 Mpc by the blue 
dotted line, the 180B02RS1 ReBEL model with j3 = 0.2 
and r s = 1 /i _1 Mpc by the green dot-dashed line and the 
180B1RS1 ReBEL model with /3 = 1.0, r s = 1 /i _1 Mpc 
by the red dashed line. Also, we indicate the Nyquist 
scale R ~ 1.4 /i _1 Mpc again by means of the vertical 
dashed line. 



TABLE IV: Values of the deviations AS 3 and AS 4 of the 
skewness and kurtosis, at a scale of R ~ 1.4 /i~ 1 Mpc, 
measured for the ReBEL models from those of the canon- 
ical ACDM model. For the definition of AS3 and AS4 
see Eqn. (|27|) . First column: redshift z. Three additional 
columns: values for S3 (top table) and S4 (bottom table) 
for three 180 h~ Mpc ReBEL ensemble simulations, 180B- 
05RS1, 180B02RS1 and 180B1RS1. 



AS 3 



z 180B-05RS1 180B02RS1 180B1RS1 



0.0 


0.126 


-0.066 


-0.247 


0.5 


0.295 


-0.038 


-0.245 


1.0 


0.270 


-0.103 


-0.328 


2.0 


0.220 


-0.100 


-0.305 


5.0 


0.020 


-0.014 


-0.120 



A5 4 



z 180B-05RS1 180B02RS1 180B1RS1 



0.0 


0.216 


-0.154 


-0.484 


0.5 


0.698 


-0.075 


-0.471 


1.0 


0.560 


-0.219 


-0.603 


2.0 


1.050 


-0.242 


-0.550 


5.0 


0.600 


-0.045 


-0.288 



1. Deviation Scale 

Earlier, we had noted that at z = the ratio of 
the hierarchical amplitudes 5„(i?)/5 ACDM in the vari- 
ous ReBEL models to that in the ACDM model is close 
to unity on large scales, scales considerably in excess of 
the ReBEL scale parameter r s and in the order of the 
scale of transition between linear and nonlinear evolu- 
tion. When assessing this ratio for skewness and kurtosis 
at other redshifts, we notice the same trend. 

Interestingly, there is a slight but seemingly systematic 
shift in the scale at which the skewness and kurtosis ra- 
tios start to deviate significantly from unity. We observe 
that this scale gradually shifts towards larger scale as 
the evolution proceeds. When looking at the scale R 5 % 
at which the skewness of the ReBEL models differs more 
than ~ 5% from the ACDM skewness, in the case of the 
P = 1.0 ReBEL model we find that at z = 5 /j -1 Mpc it 
is only R ~ 6 ft, -1 Mpc while at z = 2 it has increased 
to R ~ 10 /i _1 Mpc (see fig. [T5|) . The observed trend is 
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FIG. 11: (color on-line)Left column: redshift evolution of the skewness ratios Sf eSEi /5^ CDM for different ReBEL models. The 
lines mark the ensembles 180LCDM (solid line), 180B-05RS1 (dashed line), 180B02RS1 (dotted line) and 180BlRSl(dotted- 
dashed line). Panels shows redshifts for z = (the top panel) to z — 5 (the bottom panel). Vertical dashed line marks the 
Nyquist scale ~ 1.4 ft -1 Mpc. Righthand panel: identical set of redshift panels for the kurtosis ratio S^ eBEL /S4 DM . The 
error-bars correspond to the la scatter of the 180LCDM ensemble. 



directly linked to the scales on which the density field 
reaches non-linearity: the hierarchical amplitudes can 
only start to deviate from the canonical ACDM values 
through the related strong mode couplings. The other 
ReBEL models display similar evolutionary trends, al- 
though the details may differ somewhat. 

At more recent redshifts, in all ReBEL models the 
growth of the deviations slows down, and at z = the 
scale is still R 5 % ~ 10 /i _1 Mpc. Despite the growing 
amplitude of fluctuations at small nonlinear scales and 
the corresponding deviations of the ReBEL moments at 
these scales, the dynamical screening mechanism does 
not lead to the spread of these deviations to scales larger 
than ~ 10 /i _1 Mpc. We may expect this, since the dy- 
namical impact of the additional ReBEL scalar force will 
be rendered insignificant for Fourier modes smaller than 



the comoving Fourier mode k s = 27r/r a [l7j. The required 
strong mode coupling will therefore not materialize. This 
observation is in agreement with the behaviour of the 
power spectrum P(k) of the density perturbations, as 
noted in [l8|,[20|]- Figure [LH illustrates the convergence of 
the deviation scale R$% in the case of all three ReBEL 
models. 



2. Redshift Dependence 

Another interesting question with respect to the devi- 
ations of the ReBEL model skewness and kurtosis from 
the ACDM models concerns the issue at which redshift 
these are expected to be optimal. To address this issue, 
we assess the ReBEL S3 and S4 deviations on a scale of 
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FIG. 12: (color on-line) Skewness ratio sg eBEL / S£ CJ3M (lefthand panel) and kurtosis ratio Sf eBEL /S^ CDM (righthand panel) 
as a function of redshift z. Plotted are the ratios for three different ReBEL models: 180LCDM (solid line), 180B-05RS1 (dashed 
line), 180B02RS1 (dotted line) and 180BlRSl(dotted-dashed line). The error-bars correspond to the lcr scatter of the 180LCDM 
ensemble. 



~ 1.4 /i _1 Mpc scale. At this scale the find the highest 
deviations within the range set by the Nyquist scale. 

In table IIVI we list the values of the skewness and 
kurtosis deviations AS3 and AS4, determined for the 
180 h- x Mpc ReBEL simulation ensembles 180B-05RS1, 
180B02RS1 and 180B1RS1. The magnitude of the hier- 
archy amplitude deviations AS n is defined as: 



A5„ 



eBEL 



cACDM 



1 



(27) 



Interestingly, the most pronounced discrepancies be- 
tween the ReBEL and the standard model DM skewness 
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FIG. 13: (color on-line)Evolution of skewness deviation scale 
R s % as a function of redshift. At the scale R$% the skew- 
ness for the ReBEL model deviates by 5% from the value 
for the canonical ACDM model. Plotted are the deviation 
scales for three different ReBEL models: 180LCDM (solid 
line), 180B-05RS1 (dashed line), 180B02RS1 (dotted line) and 
180BlRSl(dotted-dashed line). 



and kurtosis are not found at the present epoch z = 0. 
Instead, we find the maximal deviations in the range 
0.5 < z < 2.0. This is directly confirmed by the visual in- 
spection of the two panels in figure 1121 where we plotted 
1 + AS 3 and I + AS4, ie. the ratios eB£t /5 3 ACDM and 
S ReBEL j S kcvM^ vcrsus re dshift z. Amongst the rather 
limited redshift archive at our disposal, the maximum 
appears to be found at z = 1. For a more precise de- 
termination of this epoch, we would need a considerably 
more densely binned redshift archive. Nonetheless, tak- 
ing into account the errors in the amplitude estimates, 
we may confidently locate the maximum somewhere in 
the range quoted at the beginning of this paragraph. 

In our numerical experiments, at z = 1 the AS3 
reaches 10.3% for the 180B02RS1 ensemble and 32.8% 
for the 180B1RS1 ensemble. The 54 deviations are con- 
siderably larger, and attain values of 21.9% and 60.3% 
for the same ensembles. We should emphasize that even 
while the deviations appear to reach their maximum at 
around z = 1, they are still substantial at the current 
epoch, attaining values of 24.7% and 48.4% for the skew- 
ness and the kurtosis in the 180B1RS1 ReBEL ensemble. 

We therefore reach the important conclusion that the 
sharpest "fingerprint" of the scalar dark matter interac- 
tions, in terms of skewness and kurtosis, should be found 
at moderately intermediate redshifts. 

The answer to the question if these signatures can be 
detected in the observational reality depends on a range 
of issues. One of the most important ones is that of the 
bias between the dark matter distribution and the baryon 
density and galaxy distribution. As yet, there is not a 
definitive insight into how much this will be influenced 
by ReBEL dark matter forces. Keselman, Nusser and 
Peebles [2(| recently studied the growth of cosmic struc- 
ture in a simulation containing dark matter and baryons. 
While they confirm the expectation that the effect of 
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RcBEL forces on the small scale baryon distribution is 
much smaller than that on dark matter, they also find 
that cannot be ignored. Preliminary results of our high- 
resolution joint simulation of baryons and DM within the 
ReBEL model (work in preparation) shows that the im- 
pact of the scalar interactions is not only imprinted on 
the moments of the dark matter density field, but also 
on the baryon density field. The implication may be that 
it will indeed be feasible to put observational constraints 
on the ReBEL cosmology parameter space with the help 
of galaxy catalogues. 

VII. DISCUSSION 

In this paper we have addressed the question in how 
far the differences between cosmological models involv- 
ing a scalar long-range dark matter interaction would 
distinguish themselves from the canonical ACDM mod- 
els in terms of their statistical properties. In answering 
this question, we focussed on the hierarchy of correlation 
functions and moments of the cosmological density fields. 

On the basis of a large ensemble of cosmological N- 
body simulations within the context of the standard 
ACDM cosmology and a range of different ReBEL long- 
range interaction models, we have attempted to identify 
the statistical differences between the models and the 
parameters and circumstances which will optimize our 
ability to discriminate between the different cosmologies. 
The simulations in this study are restricted to the pure 
dark matter distribution. 

To measure the moment and correlation functions, we 
base ourselves on the cumulants of the counts in cells 
of the particle distributions produced by the various N- 
body simulations. In the first stage of our study, we 
have thoroughly checked the accuracy and reliability of 
our implementation of the Counts in Cells method. We 
also assessed the influence of practical limitations, such 
as that of the finite length of a simulation box, on the 
measurements of the moments. This is particularly cru- 
cial given the circumstance that the gravitational force 
in the RcBEL models includes an intrinsic scale length. 
The conclusion from our experiments is that our imple- 
mentation of the CIC method succesfully recovers the 
known results from perturbation theory in the context of 
the canonical ACDM model. 

Subsequently, we have applied our toolbox to the mea- 
surement of higher order N-point correlation functions, 
from TV = 2 to TV = 9, and the related hierarchical am- 
plitudes, from the skewness S3 and kurtosis $4 to 5*8. 
Amongst the most outstanding conclusions are: 

• At scales comparable to the screening length pa- 
rameter r s of the ReBEL model, the N-point func- 
tions £at and the hierarchical amplitudes S n show 
deviations from the values expected in the standard 
ACDM cosmology. 

• In general, the amplitudes Sn become smaller as 



the DM force strength parameter (3 is larger. In 
the hypothetical situation of a negative /3, the S n 
are larger than in the case of the ACDM cosmology. 
Usually, the magnitude of the (3 = 0.2 model values 
S n are still in the order of the ACDM values, which 
technically correspond to the /3 = values. 

• In a detailed comparison between the skewness and 
kurtosis of dark matter density fields, we find that 
the relative deviation of the kurtosis in the RcBEL 
models from that in the ACDM model is consider- 
ably larger than that for the skewness. In general, 
this is true for the whole range of hierarchical am- 
plitudes: the deviation of S n in the ReBEL models 
is larger when it concerns a higher order n. 

• The deviations of the hierarchical amplitudes 
S n (R) in the RcBEL models from that in the 
ACDM model are larger at smaller scales R. At 
scales where the evolution of the density field is 
still in the linear or quasi-linear regime, the devia- 
tions are negligible. Only at highly nonlinear scales 
we notice substantial and measurable differences. 

• The scale at which we find substantial differences 
between S n (R) in the RcBEL models and in the 
canonical ACDM model gradually grows in time, 
a direct manifestation of the increasing scale of 
non-linearity as a result of cosmological structure 
growth. However, this increase comes to a halt 
when nonlinear structure growth has proceeded to- 
wards scales where the intrinsic screening length of 
the ReBEL forces calls a halt to its impact on the 
dark matter distribution. 

• The deviations of the skewness Ss(z) and kurtosis 
Si{z) of the ReBEL model from those in the ACDM 
cosmology reach their maximum in the moderate 
redshift range 0.5 < z < 2. In other words, the 
imprint of RcBEL forces in the N-point correlation 
functions should be expected to be more prominent 
at medium redshift than at the current epoch. 

By confirming that there are noticeable differences in the 
higher-order clustering patterns between the standard 
ACDM cosmology and that in the ReBEL long-range 
dark matter interaction cosmologies we have identified a 
viable path towards constraining or falsifying these mod- 
els on the basis of the observed galaxy distribution at 
moderate rcdshifts. Nonetheless, to be able to substanti- 
ate these claims we need to extend this analysis to more 
elaborate models. First, we need to assess whether the 
same significant conclusions may be drawn when the den- 
sity field is sampled on the basis of dark matter halos and 
galaxies. This is a particularly important issue as the 
small measured differences between the LCDM and the 
ReBEL models might be washed out in the observation- 
ally relevant situation where the estimates are inferred 
from the dark matter halo distribution. Also, we need to 
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investigate the extent to which these findings are influ- 
enced by working in redshift space instead of in regular 
(comoving) physical space. We foresee a substantial im- 
pact of the short-range ReBEL forces. 

In our upcoming study, we will address these questions 
on the basis of mock galaxy survey models, which will 
allow a direct comparison with circumstances prevailing 
in the observational reality. 
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